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Abstract The present work addresses the question how sam¬ 
pling algorithms for commonly applied copula models can be 
adapted to account for quasi-random numbers. Besides sam¬ 
pling methods such as the conditional distribution method 
(based on a one-to-one transformation), it is also shown that 
typically faster sampling methods (based on stochastic rep¬ 
resentations) can be used to improve upon classical Monte 
Carlo methods when pseudo-random number generators are 
replaced by quasi-random number generators. This opens 
the door to quasi-random numbers for models well beyond 
independent margins or the multivariate normal distribution. 
Detailed examples (in the context of finance and insurance), 
illustrations and simulations are given and software has been 
developed and provided in the R packages copula and qrng. 
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1 Introduction 

In many applications, in particular in finance and insurance, 
the quantities of interest can be written as where 

X = : £2 —> is a random vector with dis¬ 

tribution function H on a probability space {£2,^,V) and 
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ff): —>■ K is a measurable function. The components of X 
are typically dependent. To account for this dependence, the 
distribution of A can be modeled by 

H{x)=C{Fiixi),...,Fd{xd)), xeR‘‘, (1) 


where Fj{x) = ^ x), j G {1,.. .,£/}, are the marginal 

distribution functions of // and C: [0, l]'^ —>■ [0,1] is a copula, 
i.e., a distribution function with standard uniform univariate 


margins; for an introduction to copulas, see (McNeil et al 


2005 Chapter 5), |Nelsen ( 2006| l or Joe ( 2014) l. A copula 


model such as ([T]) allows one to separate the dependence 
structure from the marginal distributions. This is especially 
useful in the context of model building and sampling in the 
case where E['ff)(A)] mainly depends on the dependence 
between the components of X, so on C; for examples of this 
type, see Section]^ 

In terms of copula model Q, we may then write 


E['f'o(A)] = E['f'(I/)] 

where U = {U\ ,..., Ud) : 12 —is a random vector with 
distribution function C, : [0, l]'^ —> K is defined as 


'F{uu...,Ud) = %{F^ {ui),...,F^ (ud)), 

and F^{p) = inf{.x G K : Fj{x) > p}, 7 G {1, ...,£/}, are the 
marginal quantile functions. If C and the margins Fj, j G 
{1 ,...,(!}, are known, we can use Monte Carlo simulation 
to estimate E['f'(I/)]. For a (pseudo-)random sample {!/, ; 
/=!,...,«} from C, the (classical) Monte Carlo estimator 
ofE['f'(t/)] is given by 

if'f'(I/,0«E['F(I/)]. 

The main challenge of a Monte Carlo simulation is thus the 
sampling of the copula. This challenge also holds for quasi- 
Monte Carlo (QMC) methods, and is actually amplified by 
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the fact that these methods are more sensitive to certain prop¬ 
erties of the function f'. Thus the choice of the construction 
method of a stochastic representation for C can have complex 
effects on the performance of QMC methods, a feature that 
does not show up when using Monte Carlo methods. The 
present work includes a careful analysis of these effects, as 
they must be thoroughly understood in order to successfully 
replace pseudorandom numbers by quasi-random numbers 
into different copula sampling algorithms. 

Let us briefly summarize the idea behind QMC methods 
and how they can be used for copula models; more precise 
definitions on some of the concepts used here will be given 
later. The idea is to start with a so-called low-discrepancy 
point I'ef = {vi,..., v„} C [0,1)^, with k>d, which is de¬ 
signed so that its empirical distribution over [0,1)^ is closer 
(in a sense to be defined later) to the uniform distribution 
U[0,1)^ than a set of independent and identically (i.i.d.) ran¬ 
dom points is. Assuming that for U' ^ U[0,1]^ we have a 
transformation (j)c such that ^c{U') ~ C, we can then con¬ 
struct the approximation 

T'(0c(v,))«E[T^(t/)]. (2) 

Figure[^shows the points (v;) obtained using either pseudo¬ 
random or quasi-random numbers, for a transformation 0c 
designed for the Clayton copula. 

QMC methods are typically used to approximate integrals 
of functions over the unit cube via 

^ E/(v,) « / /(v)iiv = /(/). (3) 


A widely used upper bound for the integration error |/(/) — 


Qn\ is the Koksma-Hlawka inequality (see, e.g., Niederre- 


Iter 


(1992 1 ), which is of the form V {f)D* {P„), where V (/) 


measures the variation of / in the sense of Hardy and Krause, 
while D*{P„) measures the discrepancy of P„, i.e., how farP„ 
is from U[0,1)^. 

To analyze the properties of the QMC approximation Q 
for E['f'(f/)], there are two possible approaches. The first 
one is to define / = f' o 0c and work within the traditional 
framework given by ([^, the Koksma-Hlawka inequality with 
this composed function and the low-discrepancy point set P„. 
The second one is to think of Q as approximating 


E['PiU)]= nu)dC{u) 
“'[ 0 , 1 )'' 


and work with generalizations of the Koskma-Hlawka in¬ 
equality that apply to measures other than the Lebesgue mea¬ 
sure; see |Hlawka and Muck| ( |1972l l or |Aistleitner and Dick| 
( |2015| l. In the latter case, we work with the function 'P and 
view the transformation 0c as one that is applied to the low- 
discrepancy point set rather than to P. That is, here we 
work with the transformed point set = {0c (vi ),..., 0c (v„)} 




U, 

Figure 1 1000 realizations of a bivariate Clayton copula with 9 = 
2 (Kendall’s tau equals 0.5), generated by a pseudo-random number 
generator (top) and by a quasi-random number generator (bottom). 


and analyze its quality via measures of discrepancy that quan¬ 
tify its distance to C rather than comparing P„ to U[0,1)^. 

QMC methods have been used in a variety of applica¬ 
tions, but so far most of the problems considered have been 
such that the stochastic models used can be formulated using 
independent random variables (e.g., a vector of dependent 
normal variates can be written as a linear transformation of 
independent normal variates). In such cases, the transfor¬ 
mation of the uniform vector v into observations from the 
desired stochastic model can be easily obtained by trans¬ 
forming each component Vj of v using the inverse transform 
method, which is deemed to work well with QMC in part 
because of its monotonicity, and also because it corresponds 
to an overall one-to-one transformation from [0,1)'^ to K'^. 
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In the more general copula setting considered in this pa¬ 
per, at first sight the so-called conditional distribution method 
(CDM) (which is the inverse of the copula-based version of 
the Rosenblatt transform) appears to be a good choice to use 
with quasi-random numbers, as it is the direct multivariate 
extension of the inverse transform mentioned in the previous 
paragraph. Namely, it is a one-to-one transformation that 
maps [0,1 Y to [0,1 Y it is monotone in each variable. A 
transformation with k = d is certainly desirable (and prefer¬ 
able to a many-to-one transformation with k > d) when used 
in conjunction with QMC methods, since these methods do 
better on problems of lower dimension. Also, Intuitively the 
monotonicity should be helpful to preserve the smoothness 
of f' (for the first approach) and the low-discrepancy of 
(for the second approach). An additional advantage of the 
CDM is that it is applicable to any copula C (and the only 
known algorithm such general). However, the involved (in¬ 
verses of the) conditional copulas are often challenging to 
evaluate which has led to other sampling algorithms being 
more frequently used. An example is the Marshall-Olkin 
algorithm for sampling Archimedean copulas, which we also 
address in this work. 

The paper is organized as follows. Section [^provides a 
short introduction to quasi-random numbers. Sectionj^shows 
how quasi-random samples from various copulas (and thus 
multivariate models with these dependence stmctures) can be 
obtained using different sampling algorithms. Detailed exam¬ 
ples are given. In Section]^ we discuss the theoretical back¬ 
ground supporting each of the two approaches mentioned 
earlier to analyze the use of low-discrepancy sequences for 
copula sampling. Numerical results are provided in Section 
Finally, Section|^includes concluding remarks and a dis¬ 
cussion of future work. Note that most results and figures 
presented in this paper (as well as additional experiments 
conducted) can be found in the R packages copula (see the 
vignette qrng) and qrng (see demo(basket_options) and 
demo(test_functions)). 


2 Quasi-random numbers 

Here we assume that a random sample {tf, : i — 
from a copula C can be generated by transforming a random 
sample {U'j: i = 1,... ,n} from U[0,1]^ with k > t/; several 
algorithms for copula models fall under this setup. Due to the 
independence of the vectors U'l, realizations of the sample 
{U'l (obtained by so-called pseudo-random 

number generators (PRNGs)) will inevitably show regions 
of [0,1]^ which are lacking points and other areas of [0,1]^ 
which contain more samples than expected by the uniform 
distribution. To reduce this problem of an inhomogeneous 
concentration of samples, quasi-random number generators 
(QRNGs) do not aim at mimicking i.i.d. samples but Instead 
at producing a homogeneous coverage of [0,1]*. 


The homogeneity of a sequence of points over the unit 
hypercube can be measured by its discrepancy, which relates 
to the error Incurred by representing the (Lebesgue-)measure 
of subsets of the unit hypercube by the fraction of points 
in these subsets. Quasi-random sequences aim at achieving 
smaller discrepancy than pseudo-random number sequences 
and are thus also called low-discrepancy sequences. The 
rest of this section reviews concepts related to QRNG that 


are used in this paper. The reader is referred to Niederreiter 


( |1992| l, |MorokofF| ( |l998| l, l |Glasserman|2004[ Chapter 5), and 
[Dick and Pillichshammer|([2010|l for more details. 


2.1 Discrepancy 


The notion of discrepancy applies to sequences of points 
X = {vi,V2,...} in the unit hypercube [0,1)^. Denote by 
Pn = {vi,..., v„} C [0,1)^ the first n points of the sequence. 
Let be the set of intervals of [0,1)^ of the form [0,z) = 
n5=i[0,Z;), where 0 < < 1; 7 = 1! • • • I Then the discrep¬ 

ancy function of Pn on an interval [0,z) is the difference 


£([0,z);P«) 


A([0,z);P„) 

n 


-A([0,z)), 


where A([0,z);P„) =#{ i; I <i<n, v,- C [0,z)} is the number 
of points from Pn that fall in [0,z) and A([0,z)) = 115=1 
the Lebesgue measure of [0,z). 

The star discrepancy D* of is dehned by 


D*{Pn)= sup |£([0,z);P„)|. 

[0.Z)6^* 


An infinite sequence X satisfying D*{P„) G OlrT^ log^n) is 
said to be a low-discrepancy sequence. 

For a function 'P : [0,1)*^ —K, we have the well-known 
Koksma-Hlawka error bound given by 


!=1 


<V{'P)D*{P„), 


(4) 


where U' ^ U[0,1]^ and V(P) denotes the variation of the 


function f' in the sense of Hardy and Krause. See Owen 


( |2005) for a detailed account of the notion of variation and its 
applicability in practice. We also refer the interested reader to 
Hartinger et al ( 2004^ and Sobol’ (19731 for results handling 


unbounded functions (and thus of unbounded variation). 


2.2 Low-discrepancy sequences 

There are two main approaches for constructing low-discre¬ 
pancy sequences: integration lattices and digital sequences. 
Only the latter are used in this paper, so our discussion will 
focus on those. 

Digital sequences contain the well-known constractions 
of |Sobor] ( |1967| i, |Faure| ( |1982[ ), and |Niederreiter| ( |l987| i, and 
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are also closely related to the sequence proposed in |Halton| 
( |1960| l. The basic building block for this construction is the 
van der Corput sequence in base b >2, dehned as 

oo 

Sb{i) = Y, ar{i)b^''^\ i G N, (5) 

r=0 


where ar{i) is the rth digit of the Z7-adic expansion of / — 1 = 
To construct a sequence of points in [0,1)^ 


from this one-dimensional sequence, one approach is the one 
proposed by Halton (1960i, which consists of choosing k 
pairwise relatively prime integers bi,...,bk and defining the 
/th point of the sequence as 


Vi = {Sbiii),---,Sbi,ii)), iGN. 

Another possibility is to fix the base b, and choose k linear 
transformations that are then applied to the digits ar{i) from 
the expansion of i — 1 before being used in Q to define a 
real number between 0 and 1. More precisely, let Mi,... ,M<. 
be (unbounded) “oo x oo” matrices with entries in and let 


^'(0 = E E'«r+1./+!«/(*> " 

r=0/=0 


(6) 


inequality are useful to understand the behaviour of approxi¬ 
mations based on quasi-random sequences, they do not pro¬ 
vide a practical way to estimate the error. To circumvent this 
problem, an approach that is often used is to randomize a 
low-discrepancy point set in such a way that its high uni¬ 
formity (or low discrepancy) is preserved, but at the same 
time unbiased estimators can be constructed (and sampled) 
from it. Another advantage of this approach is that variance 
expressions can be derived and compared with Monte Carlo 
sampling for wider classes of functions, i.e., not necessarily 
of bounded variation (see |Owen| ( |1997a| l, [Lemieux] ( |2009| l 
and the references therein). This approach gives rise to ran¬ 
domized quasi-Monte Carlo (RQMC) methods. 

To apply this approach, we need a randomization function 
r: [0, !)■* X [0,1)^ —>■ [0,1)^ with s>k such that for any fixed 
V G [0,1)^, we have that if U' ~ U[0, !]■', then r{U',v) ~ 
U[0,1]*. Hence the individual RQMC samples have the same 
properties as those from a random sample; the difference lies 
in the fact that the RQMC samples are dependent. 

An early randomization scheme originally proposed by 
|Cranley and Patterson ( 1976|l is to take 


r{U\v) = (v-f t/') mod 1. 


( 8 ) 


where nii-j is the element in the rth row and fth column of 
Mj. Here we assume for simplicity that b is prime and all 
operations in (|^ are performed in the finite field F;,. One can 
then define a sequence of points in [0,1 by taking 

v, = {sf'{i),...,S^^ii)) (7) 


as its /th point. Sobol’ was the first to propose such a construc¬ 
tion, working in base 2 and defining the matrices Mi,... ,Mb 
so that he was able to prov e that the obtai ned sequence has 
D*{P„) e (9(n-'log^n); see ||Sobol’| ( [l96^ . 

We also point out that Halton sequences can be general¬ 
ized using the same idea as in Q. That is, one can choose 
matrices Mi,...,Mb with the elements of Mj in and 
“scramble” the digits of the expansion of / — 1 before revert¬ 
ing them via 0 to produce a number between 0 and 1. A 
very simple way to achieve this is via diagonal matrices Mj, 
each containing a well-chosen element (or factor) of Z^j. In 
our numerical experiments, we use such an approach, with 
the factors provided in Faure and Lemieux (2009 1 ; see the R 
package qrng for an implementation. 


2.3 Randomized quasi-Monte Carlo 


A randomized point set is then obtained by generating a 
uniform vector U' and letting {U') = {vi,..., v„}, where 
Vi — r{U', Vi), /€{!,...,«}. Hence the same shift U' is ap¬ 
plied to all points in P„. If we let U),..., U'g be independent 
U[0, !]■* vectors, we can constmct B i.i.d. unbiased estimators 

A' = - E 

" v,eA(t/;) 

for E['f'(f/)], whose variances can be estimated via the sam¬ 
ple variance of A„*, ■ ■ ■, A;f ■ 

In addition to the simple random shift described in (|^, 
several other randomization schemes have been proposed and 
studied. A popular randomization method for digital nets is to 
“scramble” them, an idea originally proposed by |Owen| ( |1995| l 
and subsequently studied by |Owen| ( | 1997 a| l, |Owen| ( |T997b[ ), 
|Owen|p003| l, |Matousek| ( |1998| l and |Hong and Hickernell| 
( |2003| l, among others. 

A simpler randomization for digital nets is to use the dig¬ 
ital counterpart of ([^, where instead of adding two real num¬ 
bers modulo 1, we add (in Zj,) the digits of their base b ex¬ 
pansion. That is, for u = ^ '' = 

we let 


In contrast to the error rate for Monte Carlo meth¬ 

ods based on PRNGs, approximations based on QRNGs have 
the advantage of having an error in (9(n^* log^n) when the 
function of interest T' is of bounded variation. However, in 
practice it is also important to be able to estimate the corre¬ 
sponding error. While bounds such as the Koksma-Hlawka 


M0fo V = E((“r + Vr) mod f?)/? ' 

r=0 

and define r{u, v) —u(Bb v, where the 0^ operation is applied 
component-wise to the k coordinates of u and v. The same 
idea can be applied to randomize Halton sequences (as shown, 
e.g., in Lemieux l|2009|l), but where a different base b is used 
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in each of the k coordinates. Digital shifts for the Sohol’ and 
generalized Halton sequences are available in our R package 
qrng. 


3 Quasi-random copula samples 

Sampling procedures for a tZ-dimensional copula C can be 
viewed as transformations 0c : [0,1]* —^ [0, l]'^ for some k > 
d, such that, for U' - U[0, l]'^, U := (^c{U') ^ C; that is, 
0c transforms independent U[0,1] random variables to d 
dependent random variables with distribution function C. 

The case k = d mostly known and applied as condi¬ 
tional distribution method (CDM) and involves the inver¬ 
sion method for sampling univariate conditional copulas (al¬ 
though, for example, for Archimedean copulas another trans¬ 
formations 0c with k = d is known; see |Wu et al| ( |2006] l). 
This approach thus naturally uses d independent U[0,1] ran¬ 
dom variables as input. The case k> d (often: k > d) is 
typically known as stochastic representation and is usually 
based on sampling k univariate random variables from ele¬ 
mentary probability distributions, as we will see in Section 

Ell 

In what follows we consider the above two approaches 
and show how they can be adapted to quasi-random number 
generation. 


To find the conditional copulas C{uj \ ui,...,uj^i), for j G 
{2,...,d}, for a specific copula C, the following result (which 
holds under mild assumptions) is often applied. A rigorous 
proof can be found in ( |Schmitz|2003| p. 20), an implementa¬ 
tion is provided by the function rtrafoO in the R package 
copula. The corollary that follows is an immediate conse¬ 
quence of Sklar’s theorem, for example. 

Theorem 2 (Computing conditional copulas) Let C be a 

d-dimensional copula, which, for d >3, admits continuous 
partial derivatives with respect to the first d —I arguments. 
For j G {2, ...,d} and m; S [0,1], / € {1,..., j}. 


C{uj\ui,...,uj^i) 


^j-i..AC{ui,...,Uj) 

Dj-i...iC{ui,...,Uj) 

c{ui,...,Uj_i) 


(9) 


where Dj_i,,,i denotes the derivative with respect to the first 
j arguments, C{u\,... ,Uj) denotes the marginal copula 
corresponding to the first j components and c(mi ,..., M;-i ) 
denotes the density of C{u\,... If C admits a density, 

then 0 equals 


C{uj\ui,...,uj^l) 


c{ui,...,Uj^i,Zj)dzj 


( 10 ) 


3.1 Conditional distribution method and other one-to-one 
transformations {k = d) 

3.1.1 Theoretical background 

The only known general sampling approach which works for 
any copula is the CDM. For y € {2,..., t/}, let 


Corollary 1 (Conditional copulas for general multivari¬ 
ate distributions) Let H be a d-dimensional absolutely con¬ 
tinuous distribution function with margins Fi,... and 
copula C. For j G {2,... ,d} and € [0,1], 1 G {1,.. 

C{uj I Ml,... ,My_i) =H{Ff {uj)\Ff{ui),... ,Fji^{uj^i)). 

( 11 ) 


C(My|Mi,...,My_i) =P(t/y < Uj\Ul = Ui,... ,Uj-l = Uj^i) 


denote the conditional copula ofUj at uj given U\ = u\,..., 
Uj-i = My_ 1 . If (My IM 1 , ...,My_ 1 ) denotes the correspond¬ 
ing quantile function, the CDM is given as follows; see |Em-| 
brechts et al ( 2003| l or ( Hofert|2010 p. 45). 


Theorem 1 (Conditional distribution method) Let C be a 

d-dimensional copula, U' ^ U[0,1]“^, and 0^^*^ be given by 


Ui=Uf 

t/2=C-(t/^|t/l), 


Ud=C-{U'a\Uu...,Ud-x). 

Then U = (Uu...,Ud) = 0^^"(f/') - C. 


3.1.2 Examples 

We now present several copula families and show how the 
corresponding conditional copulas and their inverses can be 
computed. To the best of our knowledge, several of these 
results have not appeared in the literamre before. 


Elliptical copulas 


An elliptical copula describes the dependence structure of 
an elliptical distribution; for the latter, see Cambanis et al| 
( [19811 l, |Embrechts et al| ( |2002| l, |Embrechts eFal ( 2003| l, or 
( McNeil et al||2005 Sections 3.3, 5). The most prominent 
two families in the class of elliptical copulas are the Gauss 
and the t copulas. 
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Gauss copulas. Gauss copulas are given by 




Proposition 1 (Conditional t copulas and inverses) With 
the notation as in the Gauss case, the conditional t copula at 
Uj, given u\,..., Uj_\, and its inverse are given by 


where <Pp denotes the li-variate normal distribution func¬ 
tion with location vector 0 and scale matrix P (a correla¬ 
tion matrix) and 0^Ms the standard normal quantile func¬ 
tion. Consider the dimension to be j and let X ~ 0p with 
X = {Xi.Q_i'^,Xj). Furthermore, assume 


V PjMJ-i) ^ii ) 


to be positive dehnite. It follows from ( Fang et al|1990 p. 45 
and 78) that 




where 


Mj|l:0-1)(-*1:0-1)) --P/',1:0-1) (^l:a-l).l:0-l)) ^-^bO'-l)’ 

( 12 ) 


C(uj I Ml, , My_ 

{uj |Ml,...,My_ 


i) = fv+,/-i(ii(y^fv'(«i)+^2)), 

, , f;+;-l(“2)Al-«2' 

\)=h 


for ii, i 2 as given in the proof. 


Figure 1^ displays 1000 samples from a t copula with 
three degrees of freedom and correlation parameter p = 
Pi 2 = Ijs/T. (Kendall’s tau equals 0.5), once drawn with 
a PRNG (top) and once with a QRNG (bottom). We can 
visually conhtm in this case that the low discrepancy of the 
latter is preserved. How this seemingly good feature trans¬ 
lates into better estimators of the form (|^ will be studied 
further through the theoretical results of Section]^ and the 
numerical experiments of Section]^ 


Archimedean copulas 


so H{xj |xi,... ,Xj-i) is again normal. With 0 ' («!:(;- 1 )) = 
(mi), ..., (m,'_i)), it follows from o that 

C{uj\ui,...,uj^i)=H{0-\uj)\0-\u^j_i))) 

0 ~Mi|l:0'-l)(^ *(“l:(i-l))) 

s/PjlW^) 

and thus that 

C^(m^' |Ml,...,My_l) 

= ^(^f|l:(7-l)(^ ^(“l:(;-l))) + iy^i|l:a-l)^ '(«;'))• 

An implementation of this inverse is provided via rtraf o (, 
inverse=TRUE) in the R package copula. 

t copulas, t copulas are given by 

Cyp(^u') typify {u\f...,ty (M J ) ) , 

where fv,/> denotes the cZ-variate ty distribution function with 
location vector 0 and scale matrix P (a correlation matrix) 
and fy * is the standard ty quantile function. The following 
proposition guarantees stability of the t copula upon con¬ 
ditioning; see the appendix for its proof and rtrafoO for 
an implementation. We are not aware that this result has ap¬ 
peared before. Given the importance of t copulas in practice, 
this is rather remarkable. 



An (Archimedean) generator is a continuous, decreasing 
function xf: [0,°°] —)■ [0,1] which satishes (/(O) = 1, i/(°o) = 
limf^K, xj/f) = 0, and which is strictly decreasing on [0, inf{f; 
xj/f) = 0}]. A -dimensional copula C is called Archimedean 
if it permits the representation 

C(«) = xi/{xi/-\ui) + --- + xi/-\ud)), 


where u= (mi ,..., u^) € [0,1]'^, and for some generator xj/ 
with inverse xf^^ : [0,1] —>■ [0,°°], where v/^^(0) = inf{f : 
xj/{t) = 0}. For applications and the importance of Archime¬ 
dean copulas in the realm of hnance and insurance, see, e.g.. 


[Hofert et al| ( |20T3] l. 


McNeil and Neslehova (2009|l show that a generator de- 


hnes an Archimedean copula if and only if xj/ is d-monotone, 
meaning that xf/ is continuous on [0,oo], admits derivatives 
up to the order I = d — 2 satisfying (—> 0 
for all / e {0,...,<i-2}, f e (0,°o), and (-1)'^- 
is decreasing and convex on (0,oo). 

Assuming xj/ to be sufficiently often differentiable, condi¬ 
tional Archimedean copulas follow from Theorem[^and are 
given by 


C(My|Ml,...,M;_l) 


^K^LxV' H«/)) 


where m/ £ [0,1], I G {I,... ,j}, and thus 


(13) 


C {uj\ui,...,uj^i) 

(14) 
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copula, then = (-1)^(1 + V^)- 

Therefore, ([T3]l equals 


C{uj\ui,...,uj^i) 
and ([T^ equals 


\2-j+LC>7V 


C {uj\ui,...,Uj^i) = 

Figure [T] displays 1000 samples from a Clayton copula with 
9—2 (Kendall’s tau equals 0.5), once drawn with a PRNG 
(top) and once with a QRNG (bottom). 


u, 



0.0 0.2 0.4 0.6 0.8 1.0 

U, 


Figure 2 1000 realizations of a f copula with three degrees of free¬ 
dom and correlation parameter p = l/\/2 (Kendall’s tau equals 0.5), 
generated by a PRNG (top) and by a QRNG (bottom). 

The generator derivatives and their inverses 

can be challenging to compute. The former are known explic¬ 
itly for several Archimedean families and certain generator 
transformations; see |Hofert et al| ( [2()12[ ) for more details. To 
compute the inverses, one can use numerical root-hnding on 
[0,1]; see rtrafo(. . . , inverse=TRUE) in the R package 
copula. This can be applied, e.g., in the case of Gumbel 
copulas. 

The following example shows the case of a Clayton cop¬ 
ula family, for which ([T^ can be given explicitly and thus 
where the CDM is tractable; this explicit formula is also 
utilized by rtraf o ( , inverse=TRUE). 

Example 1 (Clayton copulas) If !//(?) = (1 -f f > 0, 

0 > 0, denotes a generator of the Archimedean Clayton 


Marshall-Olkin copulas 

A class of bivariate copulas for which C^(m 2 |mi) is ex¬ 
plicit is the class of Marshall-Olkin copulas C(mi,M 2 ) = 
minluj “‘m 2 ,miM 2 Cti, 0:2 G (0,1), where one can show 
that 

C^(u2\ui) 

J^M2, if M2 G 

M“‘/“^ if M2G ((l-ai)M“‘('/“2-'\ 

M 2 , if M2 G '\i]- 

Figure 1^ shows 1000 samples, once drawn from a PRNG 
(top) and once from a QRNG (bottom). Here again we can 
visually conhrm the low discrepancy. 

Another class of copulas not discussed in this section 
which is naturally sampled with the CDM and thus can easily 
be adapted to construct corresponding quasi-random numbers 
is the class of pair copula constructions; see, e.g., |Kurowicka| 
|and Cooke| ( |2007[ |. For this purpose, we modihed the function 
RVineSimO in the R package VineCopula (version > 1.3). 
It now allows to pass a matrix of quasi-random numbers to be 
transformed to the corresponding samples from a pair copula 
construction; see the vignette qrng in the R package copula 
for examples. Note that if sampling of the R-vine involves 
numerical root-hnding (required for certain copula families), 
the corresponding numerical inaccuracy may have an effect 
on the low discrepancy of the generated samples. 

3.2 Stochastic representations (k > d, typically k> d) 

3.2.1 Theoretical background 

As mentioned above, pair-copula constructions are one of the 
rare copula classes for which the CDM is applied in practice. 
For most other copula classes and families, faster sampling 
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Figure 3 1000 realizations of a Marshall-Olkin copula with Ui = 0.25 
and a 2 = 0.75 (Kendall’s tau equals roughly 0.23), generated by a 
PRNG (top) and by a QRNG (bottom). 


variables. A random vector U ^ Cp^ thus admits the stochas¬ 
tic representation <P(X) = <P{AZ) for Z = (t/f),..., 

^ U[0,1]'^; here <P is assumed to act on 

AZ component-wise. Note that for Gauss copulas, this sam¬ 
pling approach is equivalent to the CDM. 

t copulas. A random vector X ~ fv,/> admits the stochastic 
representation X — sfwAX where A and Z are as above and 
W —\/r forF following a Gamma distribution with shape 
and rate equal to v/2. A random vector U ^Cy p thus admits 
the stochastic representation tv{X) = fv(-\/WAZ); as before, 
ty is assumed to act on AX component-wise. Note that 
for f copulas with hnite V, this sampling approach is different 
from the CDM. 


Archimedean copulas 


The conditional independence approach behind the Marshall- 
Olkin algorithm for sampling Archimedean copulas is one 


example for transformations 0c for A: > t/; see Marshall and 


Olkin (19881. For this algorithm, k = d +l and one uses 


the fact that for an Archimedean copula C with completely 
monotone generator \j/. 


U = iwiEi/V),...,w{E,/V))^C, (15) 

where V ^ F = independent of E\,...,Ed ~ 

Exp(l); here, F = [y/] denotes the distribution func¬ 
tion corresponding to i)/ by Bernstein’s Theorem [•] 

denotes the inverse Laplace-Stieltjes transform). To give an 
explicit expression for the transformation 0c = 0^° in this 
case, we assume that vi is used to generate V via the inversion 
method, and V;+i is used to generate Ej, for y S {1,.. .jdj. 
Then we have that 0^° = (0^°) • • • > 0^°)’ where 

(16) 


algorithms derived from stochastic representations of £/ ^ C 
are known, especially for d ^2. They are mostly class- and 
family-specihc, as can be seen in the examples below. 

3.2.2 Examples 
Elliptical copulas 

Gauss and t copulas are typically sampled via their stochastic 
representations. 

Gauss copulas. A random vector X ^ ^p admits the stochas¬ 
tic representation X — AZ where A denotes the lower trian¬ 
gular matrix from the Cholesky decomposition P = AA^ 
and Z is a vector of d independent standard normal random 


We can use a low-discrepancy sequence in A: = -|- 1 di¬ 

mensions to produce a sample based on this method. Having 
k — d+l instead of k = d isa slight disadvantage, since it is 
well known that the performance of QMC methods tends to 
deteriorate with increasing dimensions. 

To explore the effect of the transformation 0c on P,„ we 
generated 1000 realizations of a three-dimensional Hal ton 
sequence; see the top of Figure]^ where we colored points 
falling in two non-overlapping regions in [0,1)^. The first 
two of the three dimensions are then mapped via 0^™ (see 
the bottom of Figure]^ to a Clayton copula with parameter 
9=2 (such that Kendall’s tau equals 0.5). As we can see, 
the non-overlapping colored regions remain non-overlapping 
after the one-to-one transformations have been applied. To 
study the effect of the Marshall-Olkin algorithm, we look at 
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when the first dimension of the Halton sequence is mapped 
to a Gamma F(l/0,1) distribution by inversion of yi(the 
distribution of V in ( [T5] l for a Clayton copula) and the last 
two to unit exponential distributions (by inversion of 1 — vj 
for 7 = 2 ,3, so that the obtained uj is increasing in each of vi 
and Vj+i for j = 1,2). The top of Figureshows the second 
and third coordinates of the Halton sequence, and colors the 
points belonging to two different three-dimensional inter¬ 
vals (this is why not all two-dimensional points are coloured 
in the two-dimensional projected regions). We see on the 
bottom of Figure]^ that here again, the colored regions re¬ 
main non-overlapping. However, it should also be clear that 
two points in a given interval defined over the second and 
third dimension could end up in very different locations after 
this transform, if the corresponding first coordinates are far 
apart. Hence, the fact that the Marshall-Olkin transform uses 
k = d+l uniforms (and thus is not one-to-one) makes it more 
challenging to understand its effect when used with quasi¬ 
random numbers. On the other hand, because it is designed so 
that the first uniform vi is very important, it may work quite 
well with QMC since these methods are known to perform 
better when a small number of variables are important (i.e., 
see [Pick and Pillichshammer| ( |2010[ ); |Lemieux1 ( |20091 l). This 
combination (QRNG with the Marshall-Olkin approach) is 
studied further in Section]^ with numerical results provided 
in Section |5] 

Marshall-Olkin copulas 

Bivariate (d — 2) Marshall-Olkin copulas C also allow for a 
stochastic representation in our framework 0c for k> d. For 
example, it is easy to check that for (t/(, t/^) U[ 0 , 1 ]^, 

(max{t/j' , 1 / 3 “' }, max{t/ 2 ^ }) ^ C. 

This construction can be generalized to d >2 (but we omit 
further details about Marshall-Olkin copulas in the remaining 
part of this paper). 



0.0 0.2 0.4 0.6 0.8 1.0 

U, 



Figure 4 1000 realizations of the first two components of a three- 
dimensional Halton sequence with colored points in the regions 
[0, \/l78]^ and [1 — \/l78,1]^ (top): corresponding (t^M.i-faj^sfonned 
points (bottom) to a Clayton copula with 6=2 (Kendall’s tan equals 
0.5). 


3.3 Words of caution 


The plots showing copula samples obtained from QRNGs 
that we have seen so far have been promising, in that the ad¬ 
ditional uniformity (or low discrepancy) compared to pseudo¬ 
sampling was visible. Here we want to add a word of cau¬ 
tion to the effect that it is crucial to work with high quality 
quasi-random numbers, as defects that exist with respect 
to their uniformity on the unit cube will translate into poor 
copula samples. Figure illustrates this by showing two- 
dimensional copula samples obtained from quasi-random 
numbers of poor quality, corresponding to the projection on 
coordinates (20,21) of the first 1000 points of the Halton 


sequence (top) and a similar sample obtained from a gen¬ 
eralized Halton sequence (bottom), which was designed to 
address defects of this type in the Halton sequence. More 
precisely, here the problem is that this particular projection 
is based on the twin prime numbers 71 and 73 for the base. 
Defects of this type are discussed further, e.g., in Morokoff 
[and Caflisch| ( |199^ . 
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Figure 5 1000 realizations of the second and third components of a 
three-dimensional Halton sequence with colored points corresponding to 
the regions [0,0.5]^ and [0.5,1]^ (top): corresponding ())^®-transfonned 
points (bottom) to a Clayton copula with 6 = 2 (Kendall’s tau equals 
0.5). 

4 Analyzing the performance of copula sampling with 
quasi-random numbers 

In this section, we discuss the two approaches outlined in the 
introduction to analyze the validity of sampling algorithms 
for copulas that are based on low-discrepancy sequences. 


Figure 6 Samples obtained from Clayton copula with 6=2 with the 
CDM method based on coordinates 20 and 21 of the Halton sequence 
(top) and the generalized Halton sequence (bottom). 

SO that we can work in the usual Koksma-Hlawka setting 
based on uniform discrepancy measures. 

Given that a copula transform (j)c = (0c.i, ■ ■ ■, 0c.d) 
regular enough, denote its Jacobian by 

I — : ■ ■ ■) ^c,d) 

d{ui,...,Ud) 


4.1 Composing the sampling method with the function of 
interest f' 

Our goal here is to assess the quality of a quasi-random 
sampling method for copula models by viewing the transfor¬ 
mation (j)c as being composed with the function ^ of interest. 


and write 

E[^{U)]= [ ^{u)c{u)du 
"'[ 0 , 1 ]'' 

= / '^"(0c(v))c(0c(v))|/0c('/’c(v))Mv. 

J [0,1] 
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In the case of (j)c = one can easily show that 

\J^ci‘l^civ))\=c{(j)civ))-\ and thus 


for each 1 < k < i < d. Furthermore, assume that for all 
I <k <l <d and {ai,..., «/} C {1, ... ,c/}, we have 


/ ^^{^c{v))dv. 

(17) 

d^'¥{ui,...,Ud) 

Uo.i] 




While the properties of the CDM approach allow one to 
directly show ([T^ in its integral form as done above, this 
equality holds more generally for any transformation (j)c : 
[0,1]^ — >■ [0, l]'^ such that (j>c{U) ~ C whenever U ~ U[0,1]^; 
see also |Caflisch| ( |1998| l; |Pillards and Cools | ( |2006l l. 

In the case where holds, one can apply the Koksma- 
Hlawka error bound Q to transformed samples. 


(19) 


Then there exists a constant (independent of n but de¬ 
pendent on'F) such that for the choice (j)c = we have 


- '!'(«,■) - E['f'(I/)] < D* (vi,..., v„)KC^‘‘'> (M‘‘/m) 


2d-l 


Proposition 2 (Koksma-Hlawka bound for a change of 
variables) Let U C, such that ( [T7| i holds, and Ui = 
for P„ = {Vi,i =l,...,n} in[0, 1]^. Then 


1 

n 




i=i 


<D*{P„)V{^o(pc). 


Note that V(T) < oo does not imply V{^ o < oo in 
general. To get further insight into the conditions required to 
have a finite bound on the integration error, we work with a 


slight variation of the above bound that is given in (Nieder 


reiter|19^ pp. 19-20) (see also ( Hlawka and Muck|1972 


(4’)) and (Hlawka 1961 (4))), where the term Vo 0 c) is 


replaced by an expression given in terms of the partial deriva¬ 
tives of "P o 0 c assuming the latter exist and are continuous. 
It is given by 




1=1 


<D*(P„)||'Po 0 c|U.i, 


where 


|9^o0clU,i = 


LL 

1=1 a 


[ 

<9'Po0c(va|, 

■ ■ ■ ji'a/j 1) 

J[o,iy 

dvai ■ ■ ■ 

dva, 


dvai---dva, 


and the second sum is taken over all nonempty subsets 
(X = {«!,...,«/} C Furthermore, the notation 1 

in "Po 0 c(vai,... ,Va,, 1) means that each variable vj with 
j {«!,...,«/} is set to 1 . 


The following proposition provides sufficient conditions 
on the functional P and on the copula C to ensure that |j Po 
0 clU,i < °° when 0 c = 0 ^™. 

Proposition 3 (Conditions to have bounded variation with 
variable change in the CDM) Assume that P has continu¬ 
ous mixed partial derivatives up to total order d and there 
exist m,M,K > 0 such that for all u G (0,1)'^, c(u) > m > 0 
and 


where Ui = (j)^^^{vi), i = 

Proof See ( |Hlawka and Muck||1972[ (11) and the remark 
thereafter). 


Remark 1 1) As we will see in the discussion preceding the 
next proposition, in general, to ensure that ||Po 0 c||;/,i < 
oo holds, a possible approach is to bound the mixed par¬ 
tial derivatives involving P and then to verify that the 
mixed partial derivatives involving 0c are integrable. As 
explained in Hlawka and Miick (19721, Condition ( [TSl ) 
ensures that the latter condition is verified in the case of 
the CDM (or Rosenblatt) transform, and avoids having 
to deal with the function 0 c and its partial derivatives. 
Unfortunately (and while It may seem easier to work with 
the conditional copulas C{uj\ui,..., uj^ i) than with 0 c), 
in many cases the copulas Involved do not have bounded 
mixed partial derivatives everywhere, with singularities 
appearing near the boundaries when one or more argu¬ 
ments are 0 or 1. A non-trivial case where we were able to 
verify ([T^ is for the Eyraud-Farlie-Gumbel-Morgenstem 


copula (see Jaworski et al (2010 1 ), assuming the param¬ 
eters are chosen so that the density c{u) and thus the 
denominator of .. .,Uj^\) is bounded away from 

0 for all MS [ 0 , 1 )''. 

2) We note that the conditions given in ( [T9l l are not the same 
as those required to prove that 


IP 


Ua = 

l=\ a 


[0.1)' 


5'P(Mai,. 

• ■ I ) 1 ) 

du(x^ ■ • 

• duaf 


t/i/oii ■ ■ -dufxi 


is bounded; in the latter case, we only need to consider 
mixed partial derivatives of order at most one in each 
variable (since the a/s are distinct). However, in ( [T9| l, the 
j3/s are not necessarily distinct. In particular, this means 
that we need to consider the partial derivative of P of 
order d with respect to each variable and make sure it is 
bounded. 


d'^Ciui\ui,...,Ui_l) 


duay ■ ■ ■ du(Xi, 


<M, ai,.G {1,.(18) 


Let us now move away from the CDM method and con¬ 
sider a general transformation 0c. In order to study ||Po 
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0 clU,b we first need to decompose mixed partial derivatives 
of the form 


2) the function satisfies |'f'(M)| < °°for all u G [ 0 , 
and 


d’{Wo^c){vai,---,Va,A) 

dvai 


(9ft Ml . . .dPdUj 


<ooforall P 


( 22 ) 


in terms of ’f' and 0c separately. To do so, we follow Hlawka 


[and Muc^ ( |1972| ), as well as ( jConstantine and Savits|199^ 
Theorem 2.1), and obtain an expression for the mixed partial 
derivative of a composition of functions via the representation 


(9'yo0c(vaic--(Va,(l) 

dvat 


i<\p\<i 


dPiui .. .dP‘iuj 


LLcrU 




s=l y.k j 




( 20 ) 


where j3 G Nq and |J31 = T,j=i l^j- Here we do not specify 
over which values of yj and kj the inner sum in the above 
expression is taken: details can be found in the proof of 
Proposition]^ But let us point out that in the product over 
j, each index appears exactly once. On the other 

hand - and as noted in item of Remark [T] above - in the 
mixed partial derivative of "P, a given variable can appear 
with order larger than 1 . 

From ( |20l l, we see that a sufficient condition to show that 
^clU.i < °° is to establish that all products of the form 

* (9lbl0c,A:.(Vai,...,V„,,l) 

dP^u,...dP^u,l\ drj^Va,...drjJva, 

( 21 ) 


are in Li. 

We note that for the MO algorithm (assuming as we did in 
that vi is used to generate V and Vy+i is used to generate 
Ej), (pc.j is a function of vi and v^+i only, for y G {1,... ,(i}. 
Hence the only partial derivatives of 0c,; are nonzero are 
those with respect to variables in { vi, Vj+i}. This observation 
is helpful to prove the following result, which shows that the 
error bound obtained when using the MO algorithm has the 
desired behavior induced by the low-discrepancy point set 
used to generate the copula samples; note that it does not 
show that "P o ^>c has bounded variation. Its proof can be 
found in the appendix. 

Proposition 4 (Error behaviour for MO for continuous 

V) Let 0^® be the transformation associated with the Mar¬ 
shall—Olkin algorithm, as given in and that V ^ F is 
continuously distributed. Let P„ = {vi,i = be the 

point set in [ 0 ,used to produce copula samples via the 
transformation 0“*^ and let Ui — 0^‘^(v,). If 

I) the point set Pn excludes the origin and there exists some 
p>l such that m\n\<i<„Viy > Ijpn; 


with j3; G {0,... and |j31 < d; 

3) and the generator ij/{-) of the Archimedean copula C is 
such that 

a) ^^'{t) -\-t\j/”{t) has at most one zero t* in ( 0 ,°°) and 
it satisfies —t*ij/'{t*) < oo; and 

b) (1 — 1 / pn) is in Off) for some constant a > 0 ,' 

then there exists a constant (independent ofn but depen¬ 
dent of'P and 0^®) such that 


-Y^'P{Ui)-E['PiU)] 


(=1 


<c(‘'Hlog«)D*(P„). 


Remark 2 We note that if E[y] < oo, as is the case for Clay¬ 
ton’s copula family. Condition |3)|b)| can be easily checked 
via Markov’s inequality. In the case of the Gumbel cop¬ 
ula, V has an a-stable distribution and it can be shown that 
Piy > x) < for x>xq and for some constant c, where 
c and xq depend both on the parameters of the distribution; 
see (Nolan 2014 Theorem 1.12). Therefore P^^(l — 1/pn) 
can be bounded by a constant time n“ in this c ase (namely by 
max{xo, (cpn)^/“}). As for Condition 3)|a) one can show 
that t* = 9 and f* = 1 are the only zeros for the Clayton and 
Gumbel copulas, respectively. 


When F is discrete, we can split the problem into sub¬ 
problems based on the value taken by V. Then, in each case, 
the bounded variation condition is much easier to verify, 
because the transformation 0c given V is essentially one¬ 
dimensional as it is mapping each vj to an exponential i 
for j G {2,..., (9 -I-1}. Its proof can be found in the appendix. 


Proposition 5 (Error behaviour for the Marshall-Olkin 
algorithm for discrete V) Let 0^*^ be the transformation as¬ 
sociated with the Marshall-Olkin algorithm, as given in GH' 
and assume C is an Archimedean copula whose distribution 
function F ofV is discrete. Let P„ = {v,-, i = 1,..., n} the 
point set in [ 0 , used to produce copula samples via the 
transformation 0^® and let Ut = 0^’^(v;). If ( |22| l holds and 

1) there exists some p > I such that the point set P„ satisfies 
maxi<K„v,'j < 1 - \/pn; 

2) there exist constants c > 0 and q G ( 0 , 1 ) such that 1 — 
F{1) < cq^ for I > 1; 

then there exists a constant (independent ofn but depen¬ 
dent ofW and 0^®) such that 


” (=1 


<cW(log«)D*(P„). 
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Remark 3 We note that the Frank, Joe and Ali-Mikhail-Haq 
copulas are such that F is discrete. The condition on the tail of 
F stated in the proposition can be shown to hold for the Frank 
and Ali-Mikhail-Haq copulas, but not for the Joe copula (the 
distribution of V in this case has a Sibuya distribution, for 
which no moments exists, i.e., it has a very fat tail). 

Let us now move on to RQMC methods. We already men¬ 
tioned that an advantage they have over their deterministic 
counterparts is that much weaker conditions are required to 
provide variance expressions for their corresponding estima¬ 
tors. The following result shows that this also holds after 
composing "F with 0 c- 

Proposition 6 (Variance expression with a change of vari- 
ahles) If Pn = {vi ,... ,v„} is a randomly digitally shifted net 
with corresponding RQMC estimator jl„ — ^ L”=i ’f^(0c(''i)) 
and //'Var('f'(I/)) < oo with U C, then we have that 

Var(M„)= Y. (23) 


4.2 Transforming the low-discrepancy samples 

As mentioned in the introduction, we can think of 0c as 
transforming the point set P,, instead of being composed 
with "F. The integration error can then be analyzed via a 
generalized version of the Koksma-Hlawka inequality such 
as the one studied in |Aistleitner and Dick| ( |2015] l, which we 
now explain. 

Similarly to the Lebesgue case we define the copula- 
discrepancy function with respect to a copula-induced mea¬ 
sure Pc on an interval B (i.e., Pc{B) = F(U G B) for U ^C) 
as 

Ec{B-Pn) = ^^^^-Pc[B). 
n 

Let ^ be the set of intervals of [0,1 )'^ of the form [a^b) = 
YY‘j=i[ai,bj),'^h&xeQ<aj<bj<\. The copula-discrepancy 
Dc of Pn is then defined as 

Z)c(F„)= sup |£c(-B;F„)|, ( 24 ) 

Be/ 


where is the dual net of the deterministic net that has 
been shifted to get P„, and f{h) is the Walsh coefficient of f 
at h, while 


('Fo0c)(/i)= Y ^ik)P{h,k), 
kezd 


Pih,k) = 

= E 


[ 0 . 1 ) 




^2m{k-(j>c{W)-h-W) 


, W~U[0,1]"'. 


Proof It is clear from Theorem in the appendix and us¬ 
ing Representation ( [T7| l that ( |2^ holds and the condition 
Var('F(f/)) < °° with U ensures it is finite. So what re¬ 
mains to be shown is the expression for the Walsh coefficient 
of the composed function "F o 0c. It is obtained as follows: 


{W7^c){h)= [ 

= / Y 


^27ii{k,(j>c(w)b)^-2ni{h.w)i,^^ 


keZ'l 

E nt) j 


/7ci{{k,(l>c{w))b-{h~w)b) 


ke/ 


[ 0 . 1 )'' 


= Y 'Pik)Pih,k), 

ke/ 


where the third equality holds thanks to Fubini’s theorem. 


By adding assumptions on the smoothness of "F and thus 
on the behavior of its Walsh coefficients, one could obtain 
improved convergence rates for the variance given in ( | 2 ^ 
compared to the 0{\/n) we get with MC, something we plan 
to study in future work. 


and similarly for Df{P„), the star-copula-discrepancy func¬ 
tion when the sup in ( |24l l is taken over instead. 

The generalization of the Koksma-Hlawka inequality 
studied in ( |Aistleitner and Dick|2015 Theorem 1) then pro¬ 
vides 


^-YW{ui)-E[W{U)] 


<V{W)D*c{uu...,Un), 


where we assume m, = 0c(v,), i G {1,... ,n}. To get some 
insight on this upper bound, we need to know how £>c(“i > ■ • • ’ 
«„) behaves as a function of n. Unfortunately, in general we 
cannot prove that D*(vi,... ,v„) G log'^n) implies that 
Df{ui,. .. ,M„) G 6 >(n^ * log'^n). Here are a few things we can 
say, though. 

First, an obvious case for which discrepancy is preserved 
is when 0 c maps rectangles to rectangles, because then 
0c(F) G ^ for all B G and thus 

Dc{Pn)<D{Pn), 

D*c{Pn)<D*{Pn), 


where F„ = {mi ,..., m„}. However, this only happens when 
C is the independence copula, and in this case the equality 
holds. This is not a very interesting case since our focus here 
is on dependence modelling. 

For the more realistic setting where 0c does not map 
rectangles to rectangles, the following result from |Hlawka| 
|andMuck| ( |1972^ holds and gives a much slower convergence 
rate for D^(F„). 


Proposition 7 Let C be such that the Rosenblatt transfonn 
0 (G' is Lipschitz continuous on [0,1]"' w.r.t. the sup-nonn 
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II • ||oo, and {Ui = 0c(v,)} for some sequence of points {v;} 
in [0,1]'^. Then 

for some function c{d), constant in n. 


Note that the above results fully depend on the properties 
of 0c. The aim would then be to choose 0c such that a low- 
discrepancy sequence {0c('’;)} w.r.t. the copula measure Pc 
results whenever applied to a low-(Lebesgue-)discrepancy 
sequence {v,}. A more fundamental approach would be to 
directly produce a low-discrepancy sequence {«,} where the 
discrepancy is measured w.r.t. the copula measure C. This is 
something we intend to study in future work. 

Now, computing Dc or is usually not feasible in prac¬ 
tice. If we replace the sup-norm by the L2-norm, we obtain 
L2-discrepancies which are usually more practical to com¬ 
pute. Let L2-discrepancies 7 c(mi ,...,M„) and (mi,..., u„) 
be defined by 


7’c(mi,...,m„) 


'{{y,z)e[ 0 ,lp‘‘-yi<Zi}\ n ) ) 


and 


T*( \ ( [ f A{[0^z)-,Pn) . 

7’c(mi,...,m„)= / ——-C(z) 


2 ^ 1/2 

dz 


respectively. Proceeding similarly to Morokoff and Caflisch 
(19941, can be computed as 


1 


n n d 


= ^ L ■,«;.,))+ / c{zfdz 


i:=l/=li=l 

2 " fi fi 


n 


(1=1-'"/t.i -'"/t.rf 


C{z)dz. 


If we consider a convex combination C{ui,... ,uf) = 
Anf=i Mj + (1 — A)min(Mi ,... ,u^), X G (0,1), of the inde¬ 
pendence copula and the upper Frechet-Hoeffding bound, 
then one can compute Tf explicitly via 


n n d ^2 

(mi,. . . ,M„) = ^ ^ P|(l - max(MA,,;,M,,;)) -f ^ 


<:=!/=!/=! 


-f 


2{\-Xf 2X{\-X)d\ 


n d 


+ 


{d+l){d + 2 ) ' nti(2/+l) 

2(1 — ■^) ^ ^ ~ 
n t I \ ■ r • (^f-l-1)! 

k=l \ii = l 12/11 ' 




d-l d 

LL-- L 

/=lri = l 


{i + iy. 


5 Numerical results 

Through typical examples from the realm of finance and in¬ 
surance and a few test functions, we now illustrate in this 
section the efficiency of QRNG in comparison to standard 
(P)RNG for copula sampling. More precisely, we compare 
Monte Carlo sampling approaches with two types of QRNGs 
based on randomized low-discrepancy sequences: The Sobol’ 
sequence and the generalized Halton sequence, both random¬ 
ized with a digital shift. Variance/error estimates are obtained 
by using B = 25 i.i.d. copies of the randomized sequence 
and comparisons are made with MC sampling based on the 
same total number of replications. Each plot includes lines 
showing , n^ * and/or convergence rates. In addition, 
on top of each plot and for each QRNG method, we provide 
the regression estimate of a such that the variance/error is in 
0{n^^). For PRNG, we only show the results with the COM 
sampling algorithm, since the choice of method does not 
affect the error or variance very much. On the other hand, for 
QRNG we show the results both with CDM and MO (when 
applicable), as this seems to sometimes make a difference. 
Understanding better why it is so and under what circum¬ 
stances a sampling algorithm perform better when used in 
conjunction with QRNG will be a subject of further research. 

While the examples given in the next section illustrates 
the use of our proposed method in typical contexts where 
they might be used, the test functions results in the section 
that follows are meant to focus on assessing the performance 
of QRNG compared to PRNG on the sole basis of gener¬ 
ating copula samples U - without including the effect of 
the marginal distributions - and also to see if the sampling 
algorithm (CDM or MO) has an effect on the performance of 
QRNG. 

Finally, we note that QRNG based on Sobol’ point sets 
is typically slightly faster than PRNG, while the generalized 
Halton sequence runs slower than PRNG. 


5.1 Examples from the realm of finance and insurance 

Consider a random vector X = (Ai,... ,2Q) modeling d risks 
in a portfolio of stocks or insurance losses. We assume 
that the y'th marginal distribution is either log-normal with 
Xj ~ LN(log(100) + p — j G {1,... ,t/}, where 

fi = 0.0001 and a = 0.2, or Pareto distributed with the same 
mean and variance as in the log-normal case. The copula 
C of X throughout this numerical study is either a Clayton 
or an exchangeable t copula with three degrees of freedom. 
To allow a comparable degree of dependence, we will use 
the same Kendall’s tau for both models. This easily trans¬ 
lates to the parameter 0 of a Clayton copula via the rela¬ 
tionship 9 = 2t(1 — t)^' and to the correlation parameter p 
of an exchangeable t copula via p = sin(;rT/2). We denote 
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S = consider the estimation of the following 

functionals ^{X): 

■ the Best-Of Call option payoff (maxX, — K)^; 

■ the Basket Call option payoff — K)^; 

■ the Value-at-Risk at level 0.99 on the aggregated risks 

VaRo.99 (5) = F5^H0.99) = inf{x S K : Fs{x) > 0.99}, 


the expected shortfall at level 0.99 on the aggregated risks 

1 n 


ESo.99 (S) = 


1-0.99 Jo.99 ^ 


\u)du'. 


the contribution of the first and middle margin to ESq 99 
of the sum under the Euler principle, see |Tasche ( 2008| l, 

E[Xi 15 > F ^-‘ (a)] and E[Xd /2 15 > ‘ (a)]. 



1e+04 2e+04 5e+04 1e+05 2e+05 


These two functionals are referred to as Allocation First 
and Allocation Middle, respectively. 

Figures |7][^|^ and|10|(as well as Figures [13] and [T4| in 
the online supplement) display selected variance estimates 
for Clayton and t copulas with Kendall’s tau parameter equal 
to 0.2 and 0.5, using either lognormal or Pareto margins, in 
dimensions d = 5,10,20 (displayed in different rows) and 
sample sizes n G {10000,15000,... ,200000}. In the Clay¬ 
ton case, the experiment uses both the MO and CDM sam¬ 
pling methods. For the f copulas, while there is a sampling 
approach based on a stochastic representation (as seen in 
Section [3.2.2| i, there is no version of the MO algorithm avail¬ 
able, so we only use the CDM method. In addition, both 
the Sobol’ and generalized Halton QRNGs are used. In all 
cases, we see that the variances associated with the Sobol’ 
and generalized Halton quasi-random sequences are smaller 
and converge faster than the Monte Carlo variance. It is not 
clearly determined whether one sampling method is perform¬ 
ing considerably better than the other. But we note that in 
some cases, such as the estimate of the Basket Call with 
T = 0.2 in £/ = 20 dimensions (Figure]^ bottom) the MO 
sampling seems to perform better than CDM. 


5.2 Test functions 

We now consider integration results on two different test 
functions. The results are shown in Figures [TT]|12[p3| and 
[T^(the latter is in the online supplement), which are based on 
a Clayton (or f) copula with T = 0.2 and T = 0.5, respectively. 
The hrst test function is given by 

'Fi(u) = 3(uj + ... + Uj)/d, 



1e-b04 2e-i-04 5e-b04 1e+05 2e+05 


n 



1e-b04 2e-i-04 5e-b04 1e+05 2e+05 


where the vector u is obtained after transforming the uniform 
points V using either the CDM transform or the MO algorithm. 
Recall that the former requires rf-dimensional points (using 
either a PRNG or a QRNG), whereas the latter requires {d + 


Figure 7 Variance estimates for the functional Basket Call with log¬ 
normal margins based on B = 25 repetitions for a Clayton copula with 
parameter such that Kendall’s tau equals 0.2 for d = 5 (top), cl = 10 
(middle) and d = 20 (bottom). 
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Reg-coeff Sobol, GHalton & CDM: -1.69 , -1.4 Reg-coeff Sobol, GHaiton & CDM: -1.13 , -1.14 



1e+04 2e+04 5e+04 1e+05 2e+05 1e+04 2e+04 5e+04 1e+05 2e+05 



n n 

Reg-coeff Sobol, GHalton & CDM: -1.66 , -1.52 Reg-coeff Sobol, GHalton & CDM: -1.1 , -1.2 



1e+04 2e+04 5e+04 1e+05 2e+05 1e+04 2e+04 5e+04 1e+05 2e+05 



n 


n 


Reg-coeff Sobol, GHalton & CDM: -1.26 , -1.4 



1 e+04 2e+04 5e+04 1e+05 2e+05 


Reg-coeff Sobol, GHalton & CDM: -0.98 , -1 



1e+04 2e+04 5e+04 1e+05 2e+05 


n 

Figure 8 Variance estimates for the functional Best-Of Call with Pareto 
margins based on fi = 25 repetitions for an exchangeable t copula with 
three degrees of freedom such that Kendall’s tau equals 0.5 for d = 5 
(top), d = iQ (middle) and d = 20 (bottom). 


n 

Figure 9 Variance estimates for the functional VaRo ,99 with lognormal 
margins for an exchangeable t copula with three degrees of freedom 
such that Kendall’s tau equals 0.2 based on B = 25 repetitions for = 5 
(top), d = 10 (middle) and d = 20 (bottom). 
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Reg-coeff Sobol, GHalton & CDM: -1.46 , -1.1 



OD 

LU 

0) 


E 


Reg-coeff Sobol, GHalton & CDM: -1,52 , -1.41 



OD 

LU 

0) 


E 


Reg-coeff Sobol, GHalton & CDM: -1,32 ,-1.19 



Figure 10 Variance estimates for the functional ES 0.99 with Pareto 
margins for a Clayton copula with parameter such that Kendall’s tau 
equals 0.5 based on fi = 25 repetitions for d = 5 (top), d — iO (middle) 
and d = 20 (bottom). 


1)-dimensional points. Note that f'l integrates exactly to 
1 with respect to the copula-induced measure, since Uj ~ 
U[0,1], 7 S {1,..., While we know the exact value of the 
integral in this case, we still compare estimators based on B 
i.i.d. copies of either MC or RQMC, but we plot the average 
absolute error rather than the estimated variance. 

The second test function is given by 

T'2(«)=^i((0““)-'(«)), 


where 


8i{v) 


d 


n 


\Avj -l\ + aj 

l+«7 


«; = h 


which is often used as a test function in the QMC litera¬ 
ture; see, e.g., Faure and Lemieux ( 2009 ] l and the references 
therein. So here we first apply the inverse of the CDM trans¬ 
form to the copula-transformed points obtained either using 
the CDM approach or MO, and then apply the li-dimensional 
function gi. While this has the effect of simply applying the 
standard test function gi to the original sample points V; in 
the case of the CDM, in the case of the MO algorithm, we 
are not falling back on the original points v,. The hope is 
that if MO does not preserve so well the low discrepancy of 
the original points, this function would be able to detect this 
problem. 

While the second test function is mostly interesting for 
Archimedean copulas, the first one can be used more gener¬ 
ally. This is why in the results reported in Figures [TT| and \i2\ 
we also consider an exchangeable t copula with three degrees 
of freedom and Kendall’s T equal to 0.2. 

For both test functions, we see that the SoboT and gen¬ 
eralized Halton sequences always clearly outperform Monte 
Carlo, with an error often converging in rather than 

the corresponding to Monte Carlo. For the first 

function "Fi, both the CDM and MO methods perform about 
the same. We believe this might be due to the simplicity of 
Fi—a sum of univariate powers of the u/s —and the fact 
that both methods perform equally well in the univariate case 
when combined with RQMC. Looking at the results for 
we see that with RQMC the CDM method performs better 
than MO, as there is no copula transformation performed in 
the case of CDM. However, RQMC with MO is still better 
than Monte Carlo, which suggests that the MO algorithm 
manages to preserve the low discrepancy of the original point 
set. 


6 Conclusion and discussion 

The main goal of this paper was to show how copula samples 
can be generated using quasi-random numbers. This is of 
interest when replacing PRNGs by QRNGs in applications 
involving dependent samples, possibly in higher dimensions. 
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Reg-coeff Sobol, GHalton & CDM: -0.98 , -0.99 




1e+04 2e+04 5e+04 1e+05 2e+05 


n 


n 




Reg-coeff Sobol, QHalfon & CDM: -1 




1e+04 2e+04 5e+04 1e+05 2e+05 


n 

Figure 11 Average absolute errors for the test functions "Fi («) = 3{«j + 
. .. + Uj)/d (top) and («)) (bottom) for a Clayton 

copula with parameter such that Kendall’s tau equals 0.2 based on 
B = 25 repetitions for d = 5; the middle plot shows results for f'l («) 
and an exchangeable t copula with 3 degrees of freedom and Kendall’s 
tau of 0.2. 


n 

Figure 12 Average absolute errors for the test functions 'Fi (u) = 3(uj + 
... + Uj)/d (top), "ft)!/) = Kc(C{u)) + 1/2 (middle), and '¥},(u) = 
gl ({(()*“°'^)“* («)) (bottom) for a Clayton copula with parameter such 
that Kendall’s tau equals 0.2 based on 5 = 25 repetitions for d = 15; 
the middle plot shows results for •f'l (u) and an exchangeable t copula 
with 3 degrees of freedom and Kendall’s tau of 0.2. 
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We have seen that different sampling approaches can be used, 
with a focus on the COM approach and, additionally for Ar¬ 
chimedean copulas, on the Marshall-Olkin algorithm. We 
have studied the error behaviour for both methods and have 
seen that in order to prove that the composed function 'f' o 
is smooth enough to satisfy the Koksma-Hlawka bound for 
the error, sufficient conditions on the function f' are that it 
must have smooth higher order mixed partial derivatives. For 
the Marshall-Olkin algorithm, we have shown that for several 
Archimedean copula families, the corresponding transforma¬ 
tion (pc is smooth enough to guarantee the good behaviour of 
the error bound. The superiority of QRNG over PRNG for 
copula sampling was illustrated on several examples, includ¬ 
ing a simulation addressing an application in the realm of 
finance and insurance. Most of the results in this paper are 
reproducible using the R packages copula and qrng. 

Some ideas for future work would be to follow-up on 
Proposition]^ and to analyze the speed of decay of the Walsh 
coefficients of the composed function 'P o (pc, based on as¬ 
sumptions on the speed of decay of the Walsh coefficients of 
•P and the properties of the sampling method (pc- 

Concerning the copula-induced discrepancy studied in 
Section |4~2] a possible avenue for future research would be to 
construct point sets that directly minimize this discrepancy, 
without first transforming a (uniform-based) low-discrepancy 
sample. In addition, proving error bounds based on the L 2 - 
discrepancy would be useful, as this discrepancy measure is 
easier to compute. Finally, numerically computing the copula- 
induced discrepancy for samples transformed either using 
the CDM or the MO algorithm would give us some insight 
as to how conservative the bound given in Proposition |7| is. 
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Appendix 

For all the randomization schemes mentioned in Section]^ in addition 
to having a simple method to estimate the variance of the corresponding 
RQMC estimator, results giving exact expressions for this variance are 
also known and typically rely on using a well-chosen series expansion 
of the function *P of interest. The following result recalls this variance 
expression in the case of randomly digitally shifted net; see |Lemieux| 
|2009) for a detailed proof. This result is used in the proof of Proposition 
I6lin Section l4.ll 

Theorem 3 (Variance for randomly digitally shifted nets) Let P„ — 

{vi,..., v„} he a randomly digitally shifted net in base b with corre¬ 
sponding RQMC estimator Ji„ given by 

" 1=1 

and assume 'VaT('P{U)) < °ofor U f/[0,1)“^. Then we have that 


where 'R{h) is the Walsh coefficient of'R at h, given by 

V{h) = / /(M)e-2’^'(''“>'’dM 

“'[ 0 , 1 )'' 

where {h,u)i, = lY.‘j=iY,T=o^jJ‘'‘iJ+l withhjj andujj obtainedfrom 
the base b expansion of hj and Uj, respectively, and fff = {h ElP '■ 
{h,Vi)i, ^Z,,\/i= l,...,n} is the dual net of the deterministic net P„ — 
{Vi,i= that has been shifted to get P„. 


Proofs 


Proof (Proo/o/Propoifrionul AssumeP = P ' p 

r^-l ^-1 ^ 

to be positive definite matrices. Corol- 


p-l _ 

~ \ P-' p-1 

laryfHimplies that 


C{Uj\Ul,...,Uj_l) — iv V ) ’ 


where 

hvp{xi,---,Xj) = 




(25) 


r(|)(v;r)i/2v1^V 
is the density function of ty p. Using the block matrix equality 

^l:(;-l),l:(7-l) “-^1:0-1), 7 (^ 7 ,7 ) ^7,1:(7-1) ^ ’ 

we have that 
jc^P 'jc 

= -*n{7-l)^l:0'-l).l:(7-l)''^l:{7-l)+^7^7,7 +^^7^1:(7-1)'^1:(7-1),7 
= -*n(7-l)(^i:{7-l).i:{7-l)) ’-^1:(7-1) +^7^7 j' +^-^J^1:(7-1)^iT(7-1).7 

+ -*'i:(7-1)^1:(7-1),7(^7.7 ) ^7.1:{7-1)'*'47-1) 

= g + k{xjf, 

where 

s = '*^M7-l) (A:(7-1),1:(7-1)) > 

k{xjf = (xjyjp^-i-S2'j , 


Ne can thus rewrite 

ly,p(Xl,...,Xj) = <3(1-1- 


hv+j-i{l{xj)), 


vhere /zy+j-i is the density of ?v+;-i and 

r{Y±fl)y/{vTJ^T)n 


, l{Xj)^k{xj)si, Si = (/ ^ 


r{i){v7t)mV\P\ 

Jsing a change of variable argument, we compute 

f^J ( v + 7— 

ljrv.p{xi,...,xj^i, zj )dzj = a \Pf} j t 


V+7-U'W 


(/(X 7 )). 


lonsequently. 


Var()i„)= ^ |■P(/t)|^ 




='v+7-1 (/fe)), 


\p\pj. 


















Quasi-random numbers for copula models 


21 


where the last equality can be seen as follows. Let be as in 

(T2j. Since |P| = andP,-|i.{j_i) = 

by jHorn and Johnson|2012| (0.7.3.1)), we then have 


1^1 = lA:(j-l)J:(;-l)l ■ 

Proof (Proof of Proposition^ We start by providing more details on 
the expression j20^ , which is given by: 

d''Fo(pciVa,,---,Vai,l) _ 
dva, ■ ■ ■ dva, 


E 

i<\p\<i 


Ml . . . 


I 


s 


E E 

S=1 {k,Y)eps{p.a) j=l 


d^i'Va, ...di'j'Va, 


where [ Pj, and the set Ps{P ,Ot) includes pairs (k. y) 

such that k is an i-dimensional vector k = (k\,.. .,kf) where each kj 6 
{l,...,ii}, and y is an ^/-dimensional vector y= (yi,...,yj) where 
each Yj is an /-dimensional vector whose entries are either 0 or 1, and 
Yfj=i Jj.i = 1 for / 6 {1,...,/}. Finally, the C y are c onstants, which are 
defined in detail in |Constantine and Savits|(l996^ , along with further 
information on the precise definition of p, (ft, y). 

As mentioned in Section |4.1| a sufficient condition to show that 
1*^^° ^dU.i < “ is to establish that all products of the form iD are in 
L \, which we recall is given by 

Wuyywy ;y dri-^va,...drr'va, ’ 

fori 6 {1,...,/} and (k.y) 6 ps(P.a). 

Recall also that for the MO algorithm, tpcj is a function of vi and 
V /+1 only, for / = 1,... ,d. Hence the only partial derivatives of (j>cj that 
are nonzero are those with respect to variables in {vi ,v;+i}. 

Now, since we assume that l |22|l ho lds, then it means we just need 
to show that the product found in l|21|l is in Li, under the conditions 
stated in the proposition. In turn, we first show that this holds if the 
following bounds hold for the mixed partial derivatives of (jic: 


and is thus in Li as long as j28j holds. 

Case 5; 1 6 / and j such that Jj i = 1 has kj + I E /. In this case, we 
can assume w.l.o.g. that /={!,...,/} and therefore the products in iD 
are of the form 


d^IpcAvuVr+l) i-r d(l)c.j(vi,Vj+l) 

dvidvr+i '^Q+1 

and is thus in Li as long as holds. 

The last part of the proof is to show that {26\ , 1^, and @ hold. 
First we study the partial derivatives involved in these expressions and 
find they are given by: 


d<l>c.i{vuV2) 

dvi 

dtpc.ijvun) 

dv2 

dvidv2 



logV 2 dx\ 
x\ dvi ’ 

) —, 

/ ^1V2 



where X\ = F * (vi ) and = 1 / f{xi ), where / is the pdf correspond¬ 
ing to F, which exists since we assumed F was continuous. Now, the 
partial derivatives with respect to either vi or V 2 are clearly non-negative 
for all Vi and V 2 . Hence it is easy to see that and l |28| l hold, because 
we can remove the absolute value inside the integrals and therefore, 
these integrals amount to take differences/sums of (()c,r(T) at differ¬ 
ent values over its domain, which obviously yields a finite value since 
(pc.ri'c) always takes values in [0,1]. 

As for the mixed partial derivative with respect to vi and V 2 , our 
assumption on ^t'{t) +t\j/"{t) implies we have at most one sign change 
over the domain of the integral. If there is no sign change, the argument 
used in the previous paragraph to handle \26) and P8| l can be used to 
show (27) is bounded. If there is one sign change, then we let f* be such 
that 


If/(t) +t\f/'{t) <0 for 0 < t < f* and i//(f) +t\l/'{t) > 0 for t* > t. 


d<l>C.l{vi = l,v;+i) 




7+1 


dvi^l < oo, 


(26) 


/o,i)' 


5^<^C.i(''1,V2) ^ dlj)c.j{vuVj+i) 


dvidv2 


n 

7=2 


dv 


'7+1 


dvidv 2 .. -dvi < oo, and 




d(t>cAvi,Vr+l = l) ( j-j d(l>cAvi,Vj+i) 


dvi 


dv 


7+1 


(27) 


dvj \dvi < <x> 
(28) 


for all / < d-l- 1, 

We have three cases to consider. 

Case 1:1^1. Then the product in pi|l is given by 


1 

n 


d(l>c.j{vi = l.Q+l) 
dvj+i 


where we assumed w.l.o.g. that / ={2,...,/-|-l},i = / and kj = j+l 

for 7 6 {1.i}. Since each term in the product depends on a distinct 

variable, the product is in Li if P6| l holds. 

Case 2; 1 6 / and j such that Jj i = 1 has + 1 ^ /. This case can be 
analyzed w.l.o.g. by assuming I is of the form / = {l,...,r, r-|-2,...,/ + 
1} for some r> 1. In that case, the products in (D ai e of the form 


l9(jlC.r(vi,Vr+l 

dvi 


1 ) 


/-I 


n 

7=1.7 A 


^<^C.7(i'lAj+l) 

dVj+l 


Then let q(v) be a function such that — log^(v)/F *(v) = t*. For in¬ 
stance, one can verify that for the Clayton copula, ^(v) = 

When integrating the absolute value of the mixed partial derivative 
d'^^c.iyi,V 2 )/dvjdvj, we get 




logV 2 \ 
Xl J 


dV2 


dvi 


2 

Jo dvi xj 


[v/'(- log^(vi )/xi) log^jvi)] dvi 


= -2fV(f*)j^ dvi = -2Cx/{t*)iogF '(vOl^. 


(29) 


Now, in most cases ( 1) is not bounded, and thus we cannot prove 
that iR o has bounded variation. However, from there we can still get 
the upper bound on the error given in the result, by using a technique 
initially developed by|Sobor|l|1973^ to handle improper integrals, and 
later by [Hartinger et al| P004^ to deal with unbounded integration 
problems taken w.r.t. to a measure that is not necessarily uniform (as 
studied in Section]?^. Note that to apply their approach more easily, we 
need to make a small change and assume that rather than generating V as 
(vi), we use (1 — vi), so that in our study of the variation above 

(via the integral of the absolute value of the mixed partial derivatives), 
the boundedness condition fails at vi = 0 instead of vi = 1. Following 
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the approach in [Hartinger et al| l |2004[ l (see their Equation (24)) and 
taking c = (1 / pn,0 ,... ,0), the integration error satisfies 


It'l'M-nnu)] 

" i=i 


<—'?'( 1 . l)+D*iP„)V,,v{'I'o(^c)+Ire., 

pn ‘ ^ 


where jj (‘P o denotes the variation of'Fo(j)(^ over [c, 1] and 




I Vo (j)^[v)dv — J '¥o^(^{y)dv 


M 

< — for some M > 0, 
pn 


since we assumed |ip(«)| was bounded. As for V[c,i](*f' ° <^c)> we can 
infer from the steps that led to l |29| l that it is bounded by a constant 
times logF“*{l — l/p«) < alog« + logc by assumption. Therefore 
there exists a constant such that i] ['¥o (p^;) < log«. 

Proof (Proof of Proposition^^ Let p/ be such that PiV = /) = p/, for 
/ > 1. LetP; = Y!k= \ Pk for / > 1 andPo = 0. We also let (|)^(v 2 ,... ,Vd+i) 
= <l>c(f’l-W^ 2 , ■ ■ ■ for Z > 1 (transformation (pc when vi generates 

the value / for V). Consider a given value of n and low-discrepancy point 
set Pn - If we use inversion to generate V, then we have that the subset 
P^ = {f;: f/-i < f;.i < P;} will be used to produce copula samples with 
V = 1. Let Hi = \pI \ and n; = npj. It is clear that if Z becomes too large, 
then Hi will eventually he 0. Let L{n) be the largest value of Z such that 
«/ > 0, and let pi —Hiln. Then we can write 


[ Wo(j)c(v)dv--'f^'Fo(j)c{Vi) 

“'[ 0 . 1 )''+' « 

i{n) / !■ r 

^p, / 'Po(p[,{v)dv 2 -..dvd+i -^‘Po()»c(v,) 

yio.lf n, ^ 

oo „ 

+ Y. Pi _'Po())/,(v)dv2...dvrf+i 

l=L(n)+l ■'[O'l) 

/ 'i'°<j>c{v)dV2...dVd+l-YY^(Polpc{Vi) 
“'[0,1)'' ni ^ 


L(n) 

<Ypi 

l=l 


+ Y Pi 

l=L(n)+l 
L(n) 


'[ 0 . 1 )‘ 


'¥oplj{v)dV2...dVd+l 


+ Y (Pl-pl) fPo(p‘.{v)dv 2 ...dvd+i 
■Z[0,1)'' 

L(n) 

— Y iPlM"tkl)) + +C{n,d), 

1=1 


where A(n,rZ), B{n,d), and C{n,d) are bounds such that 


/ 'Poij)l^{v)dv2...dvd+i-^Y.'Po(pc{Vi) 
■/[O,!)" ni ^ 

Y Pi f ,'P°<l>c{v)dv2...dvd+i 

LW+l ■'[0.1)" 


l=L{n)+l 
L{n) 


Y (pi-pi) / 'P°(l>c{v)dv 2 ...dvd+i 

,t1 “'[0.1)'' 


< A{n,d) 


< B{n,d) 


< C{n,d). 


First, by definition of D*(P„) we have \Hi — n/l < 2nD*{P„) and thus 
\Pl ~Pl\ < 2D*(P„). Hence we can take C{n,d) = 2E(|‘P(t/)|)P)*(P„). 
Similarly, we can show that Y.T=L{n)+i Pi — oan therefore 

take B(n,d) — E(\'F{U)\)D* (P„). The analysis of the expression to be 
bounded hy Ain,d) is more complicated. First, we note that under the 
assumption we have on 'P and its partial derivatives, we need to show 


that the product in l |21) is in L \, but where each (pc.kj is replaced by 
k- ^ni^ ^ given /. Since (j>c f,, is solely a function of Vkj +\, then it means 
that the only relevant products to consider are of the form 


y dvk 




(30) 


, / — logv^ . + i \ , 

in which each term is of the form — ip ( j-i — 1 which is non¬ 
negative for any ■ Using a similar reasoning to the one used in the 
proof of Proposition[^(to conclude that j26| l and j28) hold), it is easy 
to see that l |30| l is in Li. 

What remains to be done is to analyze the discrepancy of P/,. That is, 
here we are looking for a bound on sup^^^, |£(z;P,()|, where we recall 
that is the set of intervals of [0,1)'' of the form z = nf=i 
where 0 < zy < 1. So consider a given z 6 [0,1)''. Then E{z\PI) = 
Mz'>Pn)/'kl “ ■^(z). Let zi = (P/,z) and Z 2 = {Pi-i,z), which are both 
in [0, l)'l+l. Note thatA(zi;P„) —A{z 2 \Pn) ='4(z;Pi). By definition of 
D*{P„), it is not hard to see that 


(A(zi;P„)-A(z2;Pn)) 

n 


-PiHz) 


<2D*{P„) 


and therefore 


Mz-,p!,) 

Hi 


(fi 

Hi 


A(z) 


<2Z)*(P„)^. 

ni 


Using the fact that j— n; j < 2nD* (P„), after some further simplifica¬ 
tions we get that 


ni 


<4D*{Pn)-. 

ni 


Hence we can take A(n,d) =4D*{P„)d and then get p/A(n,d) < 
AL{n)D* (P„). To show that the overall bound for the integration error is 
of the form (logn)D*(P„) times a constant, we simply need to show that 
L(n) 6 O(logf!). But this follows from our assumptions on P„ and F, 
since by definition, L[n) is the largest integer such that 1 — F{L{n)) > 
l/pn but we also have 1 — F(L{n)) < hence 

l/pn < eg"")") L(n) log( !/<?) — log c < logp + log n 

and thus L(n) < (logfi + logp-|-logc)/log(l/g), as required. 


Online supplement 

Additional numerical results 

Here we provide a few additional results for the experimental setup 
described in Sectio n[^ namely for the finance and insurance examples 
(see Figures [r3] and [l^ and then for the test functions (see Figures fTS] 
and |161 . 






































Quasi-random numbers for copula models 
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Reg-coeff Sobol, GHalton & COM: -1.12 , -1.08 



Reg-coeff Sobol, GHaiton & CDM: -1.47 , -1.25 



n 


n 


Reg-coeff Sobol, GHalton & CDM: -1.13 , -1.01 



Reg-coeff Sobol, GHalton & CDM: -1.5 , -1.11 



n 


Reg-coeff Sobol, GHalton & CDM: -0.84 , -1.08 



Reg-coeff Sobol, QHalfon & CDM: , -f.05 



n 


Figure 14 Variance estimates for the functionals Allocation Middle 
Figure 13 Variance estimates for the functionals Allocation First for Pareto margins and for a Clayton copula with parameter such that 

lognormal margins and an exchangeable t copula with three degrees of Kendall’s tau equals 0.5 based on = 25 repetitions for d = 5 (top), 

freedom such that Kendall’s tau equals 0.2 based on B = 25 repetitions ^ _ jq (middle) and d = 20 (bottom), 

for d = 5 (top), d = 10 (middle) and d = 20 (bottom). 
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Req-coeff Sobol, GHalton & CDM: -1.1 , -1.01 
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Reg-coeff Sobol, GHalton & CDM: -1.1, -1.01 
Reg-coeff Sobol, GHalton & MO: -0.87 , -0.94 
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n 

Figure 15 Average absolute errors for the test functions ‘Fi (u) = 3{«j + 
.. . + «^)/ii (top) and («)) (bottom) for a Clayton 

copula with parameter such that Kendall’s tau equals 0.5 based on 
B = 25 repetitions for = 5: the middle plot shows results for "Fi («) and 
an exchangeable t copula with three degrees of freedom and Kendall’s 
tau of 0.2. 
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Reg-coeff Sobol, GHalton 8. CDM: -1 "-1.02 
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Figure 16 Average absolute errors for the test functions Fi (u) = 3(mj + 
... + Uj)/d (top), = gi (bottom) for a Clayton 

copula with parameter such that Kendall’s tau equals 0.5 based on fi = 
25 repetitions for d = 15: the middle plot shows results for Fi («) and 
an exchangeable t copula with three degrees of freedom and Kendall’s 
tau of 0.2. 











